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Abstract 

One of the important aspects of IsoGeometric Analysis (IGA) is its link to Computer Aided 
Design (CAD) geometry methods. Two of IGA’s major challenges are geometries made of boolean 
assemblies such as in Constructive Solid Geometry (CSC), and local refinements. In the present 
work we propose to apply the Additive Schwarz Domain Decomposition Method (ASDDM) on 
overlapping domains to actually solve Partial Differential Equations (PDEs), on assemblies as 
well as a means for local refinement. 

As a first step we study a collection of simple two domains problems, using the software 
package GeoPDEs, all of them converged rapidly, even on very distorted domains. For local 
refinement we implement a ’’Chimera” type of zooming together with ASDDM iterations, and 
here again we get very attractive results. A last effect is that the additive methods are naturally 
parallelizable and thus can extend to large 3D problems. We give some examples of computations 
on multi-patch bodies. 

This new implementation of the DD methods brings many interesting problems such as: the 
definition of boundary conditions on trimmed patches, choice of blending operators on the inter¬ 
sections , optimization of the overlap size, choice of preconditionning to guarantee convergence ( 
as we do not have a maximum principle here.) 
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1 Motivation: Using IsoGeometry and Domain Decomposition Meth¬ 
ods Together 

One of the main purposes of Constructive Solid Geometry is to provide CAD with the tools allowing 
the engineer to perform Boolean operations on the bodies he works with. There are very many 
engineering forms which are obtained as unions, intersections, subtractions etc. of basic geometrical 
figures, see Figure [l] for example. Even quite complicated bodies can often be represented as results 
of such set operations applied to some primitive geometrical forms like cubes, balls, cylinders etc. 
m • To solve Partial differential equations on domains defined by such solids , it seems natural to 
find methods that take advantage of the definitions of such solids. 

IsoGeometric Analysis (IGA) [33] is a relatively new numerical method of solving PDEs, derived 
from the FEM. The main features of IGA, which are very important in engineering and which FEM 
lacks, are: (1) exact representation of the geometry and (2) the same family of basis functions that 
are used in the construction of the geometry are used in the solution of PDEs, thus avoiding the 
cumbersome mesh generation step. 

The main concept of IGA can be described in the following way. Assume that we work in 3- 
dimensional environment and consider a cube, which constitutes the reference space. The geometrical 
domain is build as a continuous deformation of the reference cube, as shown in Figure [2j and 
represents the exact geometrical form. Given a finite space of basis functions on the reference 
domain we apply the geometrical mapping to them to generate the finite space of solutions, as the 
space spanned by their images. 

Since a single patch is limited, one would like to extend IGA to complex domainss, and based 
on CSG and consider a collection of subdomains . Domain Decomposition (DD) methods are then 
natural candidates for the solution of PDEs over complex domains. More for IGA the problem of 
local zooming becomes more involved since the refinement need to be introduced at the level of the 
parametric space, where the rigid structure of tensor product must be respected. ( In view of this 
restriction different methods of local zooming were proposed for IGA. Popular ways of doing that is 
the T-spline technique |20] or the Truncated Hierarchical B Splines (THB) one, [10]. This allows 
one to locally refine the mesh, but still the tensor product structure of the parametric space leads to 
non-local changes in the mesh [8], implying that many additional degrees of freedom are introduced 
around the areas being magnified, and increasing the computational costs of the solution. ). Domain 
Decomposition (DD) methods, can also bring a solution here by using local zooming or Chimera 
type methods introduced in Dougherty at all |9j. 

Domain Decomposition partitions the original domain into a set of overlapping( or just glued 
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Figure 1: CSG example of Boolean operations. 



Figure 2: An example of IsoGeometric mapping. 
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together) connected subdomains with smooth boundary and then solve the PDEs on each of them 
iteratively, see [1] for a detailed exposition. In most implementations the meshes of contiguous 
domains match or are linked by a mortar element. It has been studied for IGA in pj and |16j. 

In the present work we introduce the definition and properties of DD methods for non-matching 
meshes over overlapping subdomain meshes as well as applications to problems of local refinement 
(zooming) . For zooming the meshes of the newly constructed subdomains are independently defined 
and non-matching , but the geometry of the whole domain is unchanged. This way of solving the 
local zooming problem allows us to introduce new degrees of freedom only in the region under 
investigation, at the same time we use the IGA framework with all its advantages. 

2 Overview of this work 

In the next part we recall some definitions and notations related to B-splines, which are the corner¬ 
stones of all the IGA technique. We consider the refinement procedures such as knot insertion and 
degree elevation, which allow us to introduce local refinement and zooming to the IGA. 

As a model example we state the classical Poisson equation in its continuous form. We then 
proceed to the weak form of this problem, which gives us quite a natural way of solving the equation 
numerically, the Galerkin method, used in FEM and IGA. Then the main ideas of IGA are explained 
with all the computations needed to perform the actual calculations of the approximate solution. 
Imposing non homogeneous Dirichlet boundary conditions is not straight forward in IGA and some 
special techniques should be used to impose them. We describe two different methods for imposing 
essential boundary conditions. 

The next part introduces the family of Schwarz Domain Decomposition methods and gives their 
description. We choose to implement the Additive Schwarz Domain Decomposition Method (AS- 
DDM), that is equivalent to a block Jacobi iterative process since it may be easily parallelized and 
we detail its continuous and discretized weak forms. 

In the sixth part the Discretized Schwarz Domain Decomposition Method is described. The con¬ 
struction of exact and approximation projection operators for the imposition of Dirichlet boundary 
conditions is treated. 

The next part illustrates the Schwarz Algorithm in the one-dimensional case. We give the detailed 
overview of the matrix operator forms obtained from the ASDDM. 

We then give numerous numerical results. Based on the GeoPDEs software we have developed 
a code that implements the ASDDM using IGA and applied it to solve different one-, two-, and three- 
dimensional problems. In this section all the convergence results are provided. We also describe how 
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the code was parallelized and show three-dimensional examples, computed by it. 


We conclude by discussing further directions for this research and open problems. 


3 B-splines 

3.1 Introduction 

The most popular geometrical forms being used in CAGD are based on B-splines. Their definition 
provides us with very important geometrical properties. In particular, the degree of a B-spline curve 
is strictly related to the number of control points and they give us local control over the form of 
the curve. B-splines are not general enough to describe exactly even very simple shapes such as 
cylinders, balls, circles, etc. The notion of Non - Uniform Rational B-Splines (NURBS) extends B- 
splines and allows one to exactly represent a wide array of objects that cannot be exactly represented 
by polynomials, many of which are ubiquitous in engineering design. In the present work we will 
consider only B-Splines .We will now provide the exact definitions. All the details can be found in 

m- 

3.2 B-spline basis functions 

Let 5 be a set of m non-decreasing numbers, < £2 < • • • < Cm- The Ci s are called knots , 

the set E is the knot vector , and the half-open interval [£i,£i+i) is the i-th knot span. If a knot Ci 
appears k times (i.e., Ci = £j + i = • • • = ^j^), where fc > 1, is a multiple knot of multiplicity k 
and the corresponding knot span does not exist. Otherwise, if Ci appears only once, it is a simple 
knot. If the knots are equally spaced, the knot vector or the knot sequence is said to be uniform; 
otherwise, it is non-uniform. 

To define B-spline basis functions, we need one more parameter, the degree of these basis func¬ 
tions, p. The i-th B-spline basis function of degree p. written as IVj,p(£)> is defined recursively as 
follows: 


N it 0 = 


1 if &<£<&+!, 

0 otherwise. 


(3.2.1) 


AUO = -A~K N kP- 1(0 + } l+P+l } N i+ 


(3.2.2) 


Ci+p — Ci £i+p+l — Ci+1 

If the degree is zero (i.e., p = 0), these basis functions are all indicators of the the intervals they 
are defined on. 
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Figure 3: B-spline basis functions of order 0, 1 and 2 for uniform knot vector E = [0,1, 2, 3,... ]. 
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Figure 4: Quartic (p = 4) basis functions for an open, non-uniform knot vector E = 
[0,0, 0, 0, 0,1, 2, 2, 3, 3, 3,4,4,4,4, 5, 5, 5, 5, 5]. The continuity across an interior element boundary 
is a direct result of the polynomial order and the multiplicity of the corresponding knot value. 
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Figure 5: Quadratic basis functions for open, non-uniform knot vector S = {0,0, 0,1, 2, 3,4,4, 5, 5, 5} 



Figure 6: An example of a B-spline curve and its control points 

In general, basis functions of order p have p — k continuous derivatives across knot £$, where k 
is the multiplicity of the value Q as before. When the multiplicity of a knot value is exactly p, the 
basis is interpolatory at that knot. When the multiplicity is p+ 1, the basis becomes discontinuous; 
the patch boundaries are formed in this way. 

3.3 B-spline curves 

In order to define a B-spline we have to provide the following information: a set of n control points 
{.Pi}f =1 , a knot vector S = of m knots, and a degree p, such that n, m and p must satisfy 

m = n + p + 1. 

Given this information the B-spline curve of degree p defined by these control points and knot 
vector E is 

n 

C (0 = X>i,p(O p i, 

1=1 

where Ni tP (u)’ s are the B-spline basis functions of degree p. 

The point on the curve that corresponds to a knot C(£j), is referred to as a knot point. Hence, 
the knot points divide a B-spline curve into curve segments, each of which is defined on a knot span. 
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If the knot vector does not have any particular structure, the generated curve will not touch the 
first and last legs of the control polygon. This type of B-spline curves is called open B-spline curves. 
In order to make the curve interpolatory the first knot and the last knot must be of multiplicity 
p + 1. By repeating the first and the last knots p + 1 times the curve starts at Po and ends at P n . 
A knot vector of such form is called an open knot vector. Here we consider all the knot vectors to 
be open if we don’t state otherwise (e.g., it is not the case when periodic NURBS are discussed). 

In Figure [6] you can see an example of a B-spline curve. We chose a non-uniform knot vector 
S = {0,0,0,1, 2,3,4,4, 5, 5, 5} and the control points as shown in the figure. The corresponding 
parametric space and the B-spline basis functions are shown in Figure [5} 

3.4 Derivatives of B-spline 

Since the B-spline basis functions are obtained by recursion from the lower degree B-spline basis 
functions and the dependence is linear, we can obtain the recursive formulas for their derivatives in 
a similar way and represent them as a linear combination of B-spline basis functions of lower degrees, 
see (T9j. 

The derivative of a basis function is given by 

N'iAO = 7 V _c N i,P- 1(0 - 7-—W+i,p-i(0- (3.4.1) 

S i+p S i Si+p+1 Si+1 

The proof of this recurrence is obtained by induction on p. 

It is a little more complicated task to calculate the NURBS derivatives since they are fractions of 
linear combinations of the B-spline basis functions, but still the recurrence relation exists and may 
be found in the Appendix. 

3.5 Multi-dimensional B-splines 

All the definitions given above refer to curves. One can easily extend these definitions to surfaces 
and bodies by taking a tensor product of an appropriate number of knot vectors and defining the 
basis functions as the products of the corresponding one-dimensional basis functions. More exactly, 
given a control net polynomial orders p and q, and knot vectors S = {£i}^tj P+1 and 

T = {uj}™+ 9+1 , a tensor product B-spline surface is defined by 


S(£,u) = EE (3.5.1) 

i =1 j— 1 

where lVj,p(£) and Mj q (v) are the univariate B-spline basis functions corresponding to knot vectors 
E and T, respectively. 
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Figure 7: An example of a two-dimensional B-spline basis function, which is a tensor product of two 
one-dimensional B-splines 


Many of the properties of a B-spline surface are the result of its tensor product nature. The basis 
functions are point-wise nonnegative, and form a partition of unity as for any (£, v) G [£i> £ra+p+i] x 

[^1 j Vm+q+l] j 

n,m / n \ / m \ 

E N iAOM j!q (v) = (E^) E m j» = L 

*=1,7=1 Vi=l / \j=l ) 

The local support of the basis functions also follows directly from the one-dimensional functions 
that form them. The support of a given bivariate function Nij- Ptq (^,v) = Ni iP (£)Mj !q (v) is exactly 
[£i> Ci+p+l] * \Pj i V?+g+l] 

Tensor product B-spline solids are defined in analogous fashion to B-spline surfaces. Given a con¬ 
trol lattice {P i, 3 ,kYi’'jk=ii polynomial orders p, q and r, and knot vectors S = [Cl, ^ 2 , • ■ ■, Gi+p+1 ]; ^ = 
[ui, V 2 , ■ ■ ■, v m+q+ 1 ], and E = [<ti, (72,..., cr; +r+ i], a B-spline solid is defined by 

n,m,l 

S(f,w) = E iV hp(0 Af j,g(i’) L fe,r(o')Pj,i,fe- 

The properties of a B-spline solid are trivariate generalizations of those for B-spline surfaces. 

One of the most important properties of B-splines is the number of ways in which the basis may 
be enriched while leaving the underlying geometry and its parameterization unaffected. We consider 
here two different ways of doing such enrichment: knot insertion and degree elevation. We also want 
to mention that not only do we have control over the element size and the order of the basis, but we 
can also control the continuity of the basis as well. We refer to m for the details and m for the 
applications to IGA. 
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Figure 8: F-mapping and a basis function transformation 




Figure 9: Example of a F-mapping and refinement 
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3.6 The IsoGeometric Concept 


The fundamental idea behind IGA is that the basis used to exactly model the geometry will also 
serve as the basis for the solution space of the numerical method. This notion of using the same 
basis for geometry and analysis is called the Isoparametric concept, and it is quite common in 
classical finite element analysis. The fundamental difference between the new concept of Isogeometric 
Analysis and the classic concept of isoparametric finite element analysis is that, in classical FEA, 
the basis chosen to approximate the unknown solution fields is also used to define the domain . 
Isogeometric Analysis turns this idea around and selects a basis capable of exactly representing 
the known geometry and uses it as a basis for the field we wish to approximate. In a sense, we 
are reversing the isoparametric arrow such that it points from geometry toward the solution space, 
rather than vice versa. 

One of the main concepts of IGA is the geometrical mapping F. It acts on the parametric domain 

0 : 

(3.6.1) 

defining the geometrical domain in the physical space. In fact it acts on the basis functions of 
the parametric space defining both the geometry and the space of approximation functions used 
for approximating the solutions of PDEs. Figure [8] shows how the shape is constructed and the 
transformation of a basis function under the F-mapping. In the present work the basis functions 
will always be NURBS or B-splines. 

3.7 Refinement and IGA 

One of the important features of IGA is that when refining the parametric domain the geometrical 
mapping remains the same resulting in preserving the geometry while refining the solution space. 
Thus, IGA allows us to make local refinements without changing the geometrical shape. An example 
of a two-dimensional knot-insertion is shown in Figure [9] 

One of the drawbacks of IGA is the nonlocal nature of the refinement due to the tensor product 
structure of the B-splines. The difficulties arise in large scale problems since the number of degrees 
of freedom grows very fast. Different algorithms, such as T-splines [8], were introduced into IGA 
in order to avoid these problems. All of them still makes the refinement non-local since some tensor 
product structure needs to be preserved. 

This led us to think about different methods of local mesh refinement, which would allow us to 
keep the computational costs low. We propose to use the Domain Decomposition Method. Its main 
difference from the refinement techniques is that the new refined mesh is an absolutely separated 


10 


domain and does not match the coarse one, but now we have to perform an iterative procedure in 
order to obtain the solution. 

3.8 The Alternating Schwarz Method 

The Alternating Schwarz method belongs to the class of domain decomposition methods for par¬ 
tial differential equations. Domain decomposition methods can be regarded as divide and conquer 
algorithms. The main idea is the following: given an open domain (we will assume it to possess 
some additional properties in the sequel, like being smooth and connected) D we partition it into a 
number of subdomains {D*}” =1 , satisfying some conditions, and such that Q = Uf=i The original 
problem then can be reformulated as a family of subproblems of reduced size and complexity defined 
on the subdomains. 

One of the major advantages of domain decomposition methods is the natural parallelism in 
solving the subproblems defined on the subdomains. With the upcoming of parallel computer ar¬ 
chitectures, domain decomposition methods have become very popular during the last two decades. 
However the origin goes back to 1870. In [12], H. A. Schwarz introduced an algorithm to prove the 
existence of harmonic functions on irregularly shaped domains. Today this algorithm is known as 
the Alternating Schwarz Method. 

In the next section we will define the problem we are going to solve and the main tool we use to 
solve it: the IGA framework. 

4 Notations and the Problem formulation 

Given a bounded connected open domain 0 with Lipshitz continuous let us consider the Poisson 
equation 


— div(fc(x)gradu) = / 

< = A 
u — g 


in O, 
on T n , 
on v D , 


(4.0.1) 

(4.0.2) 

(4.0.3) 


where Tp (J Tw = dQ, T d and T ; y are disjoint, and n is the unit outward normal vector on <9H. The 
functions / £ 77” 1 (D),g r £ C 2 (Td) representing the Dirichlet boundary conditions and h £ C 2 (T]\r) 
representing the Neumann boundary conditions, are all given. 


For a sufficiently smooth domain, a unique solution u satisfying (4.0.1) — (4.0.3) is known to exist 
under some standard conditions. 
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4.1 Weak form of the problem 


The numerical technique we are going to take advantage of begins by defining a weak, or variational, 


counterpart of (4.0.1) — (4.0.3). To do so, we need to characterize two classes of functions. The first 


is to be composed of candidate, or trial functions. From the outset, these functions will be required 


to satisfy the Dirichlet boundary conditions of (4.0.3) 


We may now define the collection of trial functions, denoted by S as 

S = {u:u£ U 1 (Q),u\t d = g}. 

The second collection of functions in which we are interested is called the test or weighting functions. 
This collection is very similar to the trial functions, except that we have the homogeneous counterpart 
of the Dirichlet boundary conditions: 


V = {w : w E % 1 (fl),u;|r D = 0}, 


(4.1.1) 


which is a Hilbert space with respect to the associated energy norm ||it| | = a(u, u) 1 / 2 , that follows 
from the coercivity of the bilinear form a(-, •), as shown below. 

We may now obtain a variational statement of our boundary value problem by multiplying the 


equation (4.0.1) by an arbitrary test function w E V and integrating by parts, incorporating (4.0.2) 
as needed. 

Given /, g and h find u 6 S such that for all w 6 V 

/ k(x)gradw ■ gradudfl = / wfdQ + / whdT. (4-1.2) 

J J JT'n 

Define a bilinear form a(-, •) : ’H 1 (fl) x - H 1 (D) — > M and a functional L : ’H 1 (fl) — > M as: 


and 


Now the weak form reads as: 


a(w,u) = / &:(x)gradu; • gradudfl, 

Jn 


L(w) = / wfdQ + / whdT. 

J J Tjv 


a(w, u) = L(w), 


(4.1.3) 

(4.1.4) 

(4.1.5) 


We shall note here, that the bilinear functional a(-, ■) is symmetric and positive definite for any 
k(x) > 0 when considered over the space V. Indeed, as shown below it is coercive, which results in 
the symmetry and positive definiteness of the stiffness matrix we are going to define below . One 


will find in (CIARLET) that the above problem 4.1.5 has a unique solution in 'H 1 (D) 
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4.2 Galerkin’s Method 


In Galerkin’s method we define finite-dimensional approximations of the spaces S and V, denoted 
Sh and Vh such that 4 C 5, Vf, C V. These subspaces will be associated with subsets of the space 
spanned by the isogeometric basis. 

As we have just mentioned we can further characterize Sh by recognizing that if we have a given 
function (jh 6 Sh such that gh\r D = 9, then for every Uh £ Sh there exists a unique Vh £ Vh such 
that 

Uh = v h + 9h- 

This clearly will not be possible for an arbitrary function g. The problem of approximation of 
the function g by a function gh £ Sh will be treated in detail below. So let us assume at present 
that such a function gh exists. 

As in [14] we can write the Galerkin’s form of the problem as: given gh and h, find Uh = Vh + gh, 
where Vh £ 14 , such that for all Wh £ Vh 


a{wh, Uh) = L(w h ), (4.2.1) 

or 

a(w h ,v h ) = L(w h ) - a(w h ,g h ). (4.2.2) 

We will use this form for the variational formulation for the Schwarz Method. 

In general the weighting space Vh of test functions may be different than the space Vh of trial 

functions Vh , that is Vh £ Vh but Wh £ Vh / Vh, but in the following we will assume they are the 

same. 


5 General IGA 

The goal of IGA, as it is also for FEM, is the numerical approximation of the solution of partial 
differential equations (PDEs) [4]. In both approaches the PDEs are numerically solved using the 
Galerkin procedure: the equations are written in there equivalent variational formulations, and a 
solution is sought in a finite dimensional space with good approximation properties. The main 
difference between the two methodologies is that in FEM the basis functions and the computational 
geometry (i.e., the mesh) are defined element by element, using piecewise polynomial approximations, 
whereas in IGA the computational geometry and the space are defined exactly from the (input data) 
information and the basis functions (e.g., NURBS, T-splines or generalized B-splines) given by CAD. 
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Figure 10: Example of a F-mapping and the basis function transformation 


Another important difference is that the unknowns to be computed, coefficients that define the 
solution as a linear combination of the basis functions (degrees of freedom), are not nodal values, 
but “control values” . 


Let us consider a two-dimensional case where we assume that the physical domain C M 2 is 
open, bounded and piecewise smooth. We also assume that such a domain can be exactly described 
through a parametrization of the form 

F-.n->n, (5.0.3) 


where is some parametric domain (e.g., the unit square), and the value of the parameterization 
can be computed via the isogeometric mapping F, which is assumed to be smooth with piecewise 


smooth inverse, see Figure 10 


We now construct the finite-dimensional spaces Sh and Vh as spans over their basis functions: 


Sh ■= (01, 02, 4>N h ,4>N h +l, <fi Nh+N b) C S , 

and 

Vh ■ = (01, 02, 0ivj c V, 

where we first number all the basis function which vanish on the boundary r^: 0i, 02, 0jv h , and 
then the other basis functions 4>N h + 1, •••, 4>n +N b > which do not vanish on this boundary. 

The IGA prescribes us the way of building the basis functions for the spaces Sh and Vh- Let 
{4>j}j£XuB be a basis for Sh, with X U B a proper set of indices; we partition the indices in the 
same way as before: X corresponds to the basis functions that vanish on the boundary T£> and 
B corresponds to those which do not vanish on r^. With the assumption made on F, the set 
{0j o F _1 }j e i uB = {0j}j e xuB is a basis for Sh and {0j o F= {0j}jex is a basis for Vh- 
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In IGA, as introduced in [XT], the spaces Sh and Vh are formed by transformation of B-spline 
functions. We define these space in the following general way: 


S h := {u h E S : u h = u h o F E S h } = {u h € S : u h = u h o F 1 ,u h E S h }, (5.0.4) 

Vh. ■= {v h e V : v h = v h o F E Vh} = {vh e V : v h = v h o F ,v h 6 Vft}, (5.0.5) 


where t is a proper pull-back, defined from the parametrization (3.6.1), and Sh and Vh are finite 


spaces generated by the basis functions defined in the parametric domain Q. 
Hence, the discrete solution of the problem can be written as 


Uh = 


E 


atjfij — 


ctj4>j oF 1 . (5.0.6) 

j£XUt3 j£XUB 

We assume that for the function g representing the Dirichlet boundary conditions, there exists 
a function gh E Sy t such that gh\r D = 9, meaning that there exist {y j}j£B such that g = 

The more general case of g E S will be treated below. 

The solution now read as 

Uh = 7 jfij + 

j£B j&x 

Thus, instead of considering the function Uh E Sy we can treat without loss of generality the function 
Uh — J2jeB Tj'Aj ^ Vh, reducing the problem to the homogeneous Dirichlet boundary conditions case. 
So we can write 


u h = y^ j a 


jvj- 


jex 


Substituting this expression into (4.2.1), and testing against every test function <f>i E Vh, we 
obtain a linear system of equations where the coefficients aj are the unknowns, and the entries of 
the matrix and the right-hand side are a((j>i,cj>j ) and L((j>i), respectively. 


= L(4>i),Vi E 1 

iex 


(5.0.7) 

These terms have to be computed using suitable quadrature rules for numerical integration 


6 Model Problem: Poisson Equation with homogeneous Dirichlet 
BC 


We now specialize the general framework of the previous section to the particular case of the Poisson’s 


problem (4.0.1) — (4.0.3), defined in a physical domain 0 described and discretized with B-Splines. 
This constitutes the model problem on which we show in detail the IGA method. 


15 






Let us assume that the computational domain is constructed as a singleB-Spline patch, such 
that the parameterization F is given as above, see Figure |T7jj The assumptions on F of the previous 
section are here supposed to be valid. For the sake of simplicity, homogeneous Dirichlet boundary 
conditions are assumed. 


The linear equations from 5.0.7 can be rewriten as: 


N h 

y~] AijOtj = fi + hi for i = l,...,N h , 
3 =1 


where 


( 6 . 0 . 8 ) 


Aij = / A:(x)grad^grad())jdx for i,j = 1,..., N h , 


(6.0.9) 


fi= f<Pidx for i = 1,..., N h , 

J n 


( 6 . 0 . 10 ) 


hi = / hcfridr for i = 1,..., N k - 

Jr N 


( 6 . 0 . 11 ) 


Here A t j are the coefficients of the stiffness matrix, and /) and gi are the coefficients of the right- 
hand side contributions from the source and the boundary terms, respectively. All the coefficients 


are given by the values of the integrals in (6.0.8), that are numerically approximated by a suitable 
quadrature rule [la¬ 
in order to describe this rule, let us introduce K. k := {A'/ c }^ 1 , a partition of the parametric 
domain into N e non-overlapping subregions, that henceforth we will refer to as elements. The 
assumptions on the parameterization F ensure that the physical domain Q can be partitioned as 

__ N e 

n = U F (^u 

k =1 

and the corresponding elements Kk '■= F (Kk) are also non-overlapping. We will denote this partition 
by JC h :={K k }^ v 

For the sake of generality, let us assume that a quadrature rule is defined on every element /C^. 
Each of these quadrature rules is determined by a set of nk quadrature nodes: 


and by their corresponding weights 


{xi, k } C K k ,l = 1, ...,n k 


{wi )k } CR,( = 1, ...,n fc . 
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In this work we use Gaussian quadrature rules as described in TT). 

After introducing a change of variable, the integral of a generic function cj) E C^{K k ) is approxi¬ 
mated as follows 




n-k 

F(x))| det(-DF(x))|ofi? ~ ^ wi tk (j)(xi tk )\ det(AF(x ijfc ))|, 

i =1 


( 6 . 0 . 12 ) 


where xj k := F (xi k ) are the images of the quadrature nodes in the physical domain, and D F is the 
Jacobian matrix of the parameterization F. 

Using the quadrature rule, the coefficients of the stiffness matrix are numerically computed as 

-A/e 'H'k 

Aij — EE k(xi } k)wi, k grad4>j(xi tk )gr&d(pi(xi, k )\ det(DF(xi jk ))\, (6.0.13) 

k =1 1=1 


while the coefficients /, of the right-hand side are approximated as 


-/Ye 

/•-EE f (xi !k )wi,k<l>i(xi,k) I det(-DF(xj )fe ))|. (6.0.14) 

fc=l Z=1 


7 Dirichlet boundary conditions 

The imposition of homogeneous boundary conditions in IGA is straightforward. In fact, homogeneous 
Neumann conditions are satisfied naturally because they are included in the weak formulation itself. 
On the contrary, as far as non-homogeneous Dirichlet boundary conditions are concerned, a suitable 
transformation has to be applied to reduce the problem to the case of homogeneous boundary 
conditions since there is no place in the Galerkin formulation of the problem where the Dirichlet 
boundary conditions can be imposed. For this reason we must build them directly into the solution 
space. 

In FEM the imposition of Dirichlet boundary conditions can be performed easily since the basis 
functions are interpolatory and local. On the contrary, in IGA the treatment of the Dirichlet bound¬ 
ary conditions is a much more involved problem. The basis functions (B-splines or NURBS) are not 
local, nor interpolatory. 

In section [5] we assumed that there exists a function (jh E Sh such that gh\r D = 5S and refered to 
this function as a lifting. In practice, this will frequently be the case, but there will also be instances 
in which a lifting is only an approximation of g. 

The standard approach consists in constructing a function R g as suitable extension of g such 
that 

R a € n\n), 
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Rg = g on r D , 


and consider the problem (4.1.2), (4.1.5) for the difference u — R go , see [B] for details. To construct 
an approximation of R gD it is necessary to have an approximation of the boundary data function g. 
A possible strategy consists in obtaining an approximation g^ of g in the space 


Sh\Vh {^Nn+i, • • ■ j 4>i\r h+ ]srb). 


as combination of the A^-basis functions of Sh which do not vanish on To, that is 


N h +N* 

S 7 i(x) := ^2 %^|r D (x), Vx G T D . (7.0.15) 

N h +i 

Then, an approximation of R g in Sh, can be constructed as: 

N h +N b h 

R gh (x) := ^ gi0i(x), Vx G flh- (7.0.16) 


N h +1 


The control variables qi in (7.0.16) have to be computed according to the approximation method 
used . 


The easiest way is to apply these conditions directly to the control variables, as 

N h +Nl 

gh (x) := ^2 ff(Ct)^*lrx>(x), Vx G r D ; 

N h +1 

for proper values ^, i = Nh + 1,..., Nh + A r ^. Such a strategy may result in a lack of accuracy in 
case of inhomogeneous conditions and a poor accuracy of the approximated solution, even if spaces 
Vh with high approximation power are used. 

A standard approach consists in imposing Dirichlet boundary conditions by interpolation of the 
boundary data g. Another way of imposing Dirichlet boundary conditions is to use least-squares 
approximation of the given function g by a B-spline. We will review both methods. 


7.1 Least-squares approximation of the Dirichlet boundary conditions 


The least-squares approximation of the function g consists in solving an optimization problem. More 
exactly, we need to find such coefficients {Qi}2=N^+i would minimize the following integral: 


mm 

w\ Nh+Nbh 

\ l J’-n=N h +i 


N h + N h 

j (#(x) - ^2 <?i<Mx)) 2 dr. 

JVd N h +1 


(7.1.1) 
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We perform the following transformations of the integral in (7.1.1): 




( 5 '(x)) 2 dr - 2 



f Nh+N b h 

/ Y %g j 0j(x)0 i (x)dr. 
i,j=N h +1 


Let us rewrite this equation in the following form: 


where 




N h +N b h \ - 

Y ®<M X ) I dr 

i=N h +1 / 


[ ( 5 (x)) 2 dr 

Jr D 


2q 7 g + q J Mq, 


q — ( QN h +i ■ ■ ■ q Nh+N b), 


(7.1.2) 


(7.1.3) 


(7.1.4) 


/ 


g = 


J 9(*)<l>N h + 1 (x)dr 


\ 


\/5(x)^ +JV fc(x)dry 


(7.1.5) 


and M = {/ (x)<A,(x)dr}^ 2^+1 * s mass ma t 1- ix. 

In order to minimize our target function we just differentiate it once with respect to the coeffi¬ 
cients qi to get the closed form solution: 


q = M\g. 


(7.1.6) 


The matrix M being invertible since it is the Gram matrix of the basis functions of the boundary 
<9f2 which is a union of NURBS or B-spline patches in the standard £ 2 inner product. 

We consider the boundary of the physical domain T which is an image of the boundary T under 
the isogeometric mapping F and is itself a union of NURBS or B-spline patches with dimensions one 
less than the dimension of the original patch Cl. So we may decompose the boundary T = Uj 7 i into 
a union of NURBS patches and consider the minimization problem of approximating the boundary 
function on each of them. 

After we have computed the Dirichlet degrees of freedom q j, we can construct the desired function 

N h +N b h 

R 9 h( x ) := Y 9*^*( X )> Vx G Clh- (7.1.7) 

N h +1 

We now solve the homogenous problem in the variational form: 

a(w h ,u h - R gh ) = L(w h ),Vw h £ V h , (7.1.8) 

where we have only the first Nh degrees of freedom to be found from this linear system. 
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7.2 Quasi-interpolation of the Dirichlet boundary conditions 


Here we consider the point-wise approximation: we are looking for a B-spline function (on the given 
part of the boundary) which coincides with g at a number of points in the spirit of collocation. 

Consider a two-dimensional domain. Assume that we deal with a part of boundary 7 that corre¬ 
sponds to one side of the parametric domain square and there are n degrees of freedom corresponding 
to this part of boundary. We want to take some n points on the curve 7 and make the B-spline 
function we are looking for coincide with g at these point. This gives as a square system of linear 
equations: 


X/ lj4>j( x i,yi ) = 9(xi,Vi),i = 1 ...n. (7.2.1) 

3 = 1 

In order to solve this system we must require its matrix to be non-singular, and this obviously 
depends on the choice of the points (x tl ?/i)(Li • In our work we chose the n points to be the images 
of the centers of the n knot-spans of the parametric domain. It may be easily shown that the matrix 
we obtain in such a case is nonsingular. Indeed, we know that the monomials form a totally positive 
basis on the interval [0,1] and B-splines are obtained from them by the transformations saving the 
totally positivity properties [5], consequently the B-spline basis we have is also totally positive and 
the matrix we get is nonsingular. 

Thus, the system now looks like: 

n 

= 9(F(£.i,Vi)),i = 1 ■■■n, (7.2.2) 

3 = 1 

for the specific we choose. 

One of the questions we are interested in is how to choose these points optimally for the B-spline 
and NURBS geometry. 

Applying the same algorithm to all the sides of the domain where the Dirichlet boundary condi¬ 
tions are imposed we may compute all the Dirichlet degrees of freedom. 

There are also other approaches if imposing the Dirichlet boundary conditions such as local 
least-squares, see III] or Nitsche’s method [2] . 


8 Neumann boundary conditions 


Neumann boundary conditions of the form (4.0.2) are frequently referred to as natural boundary 


conditions. This is because of the way they automatically arise in the variational statement of 
a problem. Let us assume for the moment without loss of generality that k(x) = 1 and that the 
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Figure 11: Domain Decomposition 


Dirichlet conditions have already been imposed. Multiplying by a test function and integrating leads 
us to 


0 = / w(Au + f)dCl 

Jn 

= — VwVudQ + / wVu ■ ndr + / wfdVL 

Jn Jr Jn 

= — VwVudQ + / wVti • ndr + / wfd£l, 

J n J^n J ^ 


(8.0.3) 


where in the third line we have used the fact that the weighting space is defined such that w |r D =0. 
The integration by parts has naturally introduced a boundary integral over Ttv that refers ex¬ 


plicitly to the condition that we would like to impose. Using the condition (4.0.2) we simply replace 
Vu • n with the value we are imposing, h, resulting in 


— / VtcVudD + / whdT + / wfdVL = 0. 

Jr Jn 


(8.0.4) 


Now that we have formulated a discretized analog of the original problem in the next chapter we 
are going to consider a variant of Schwarz Domain Decomposition Method to solve it. 

9 Schwarz Additive Domain Decomposition Method 


9.1 Notations 

Consider the two-dimensional domain O as shown in the Figure [TTj which can be represented as a 
union of two overlapping subdomains Q = Pi U P 2 - In this work we consider only domains which 
satisfy some reasonable conditions: they are connected and piecewise smooth boundary. We want 


to solve the Poisson problem (4.0.1) — (4.0.3) on this complex domain. We denote by = 1,2 


the part of the boundary of D, which belongs to the boundary dfl of the whole domain and by 
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r \,i,j = 1,2 ,i ^ j the boundary of H*, which belongs to the second subdomain: Tl = dQi n f lj. 
The earliest known domain decomposition method is the Alternating Schwarz Method (ASM) dating 
back to 1870 [12j. The original alternating Schwarz method is a sequential method that successively 
solves the following problems: given the Dirichlet boundary conditions function g on T and starting 
with the initial guess u 2 = 0 (in particular = 0 ), iteratively obtain solutions Ui(x,y) and 

u 2 (x,y) on and ^2 respectively (n E N) [17] . 


Au" = / in Hi, 

u\ l = 9i on T 1D , 


where g 7 ?| r 2 = u 2 1 \ r 2 1 ,9i |ri = 9 Irq and 


(9.1.1) 

(9.1.2) 


A u 2 = f in O 2 , (9.1.3) 

u n 2 =g% onT 2D , (9.1.4) 

where g 2 | r i = it”| r i, g 2 |r 2 = g|r 2 - It is important to notice that the solution of the first problem 
is required before the second problem can be solved. The part of the boundary conditions of the 
second problem, g 2 = | r i, is the current solution of the first problem on the interface. It is not 
difficult to see that in this simple case of two subdomains it is an analog of a block Gauss-Seidel 
method and the two problems cannot be solved in parallel. For a straight forward parallelization 
of the algorithm an analog of the block Jacobi method is applied giving us the Schwarz Additive 
method: 


Au" = / in Hi, (9.1.5) 

Ui=9i onTi D , (9.1.6) 

where g ^\ r 2 = 5 "|rq = and 

Au 2 = f in H 2 , (9.1.7) 

u 2 =g 2 on r 2fl , (9.1.8) 

where g %| r i = u? -1 | r i, 92 |r 2 = fl'lrv 

The iterative procedure starts given two initial guesses u\ and u 2 on each subdomain. The 
boundary condition of one of the subproblems, g r - = u ™~ 1 | r ,, i. j = 1 , 2 ,i ^ j is now the solution 
of the previous iterate on the interface, hence the subproblems can be solved independently. In the 
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general case of many subdomains a great deal of parallelism can be introduced, using a method called 
multi-coloring. This is a graph theory technique that identifies each subdomain with a color such 
that disjoint subdomains have the same color; usually only a small number of colors is needed. Thus, 
the subproblems defined on subdomains of the same color can be solved independently in parallel 
and only the number of distinct colors different threads is needed to perform the computations. We 
are using such techniques below when considering 3-dimensional multi-patched domains. 

For the ASDDM as presented here, one can see that there is no notion of a global approximate 
solution. One can consider defining the global solution by choosing the subdomain iterative solution 
on each subdomain and any weighted average within the overlap: 

u n = X\v n + X2W n , 

where xi = 1 on fii\(fii 0 ^2), X2 = 1 on 0 ^2) and xi + X2 = 1 on fii 0 This approach 

extends easily to a case of multi-patched domain decomposition. 


9.2 Convergence conditions 


The iterations are performed until certain convergence conditions are met. Throughout this work we 
consider the C 2 errors as the termination condition. In the model problems when the exact solution 
is known the errors are calculated as a £? norms between the approximate solution and the exact 
one. 


Assume the exact solution u ex £ H 1 is given. Then the error is calculated as 



(9.2.1) 


We sum up the integrals on the elements and perform numerical integration element-wise. As¬ 
sume that an integration rule is given by the mesh of points {xi t k}i,k an d weights {wi^i^k- For each 
element K^ we obtain: 



(9.2.2) 
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This numerical integration formula was implemented in the code we used. An analogous formula 
for the % 1 errors can also be developed. 


10 Weak form of ASDDM 


We may now obtain a variational form of the ASDDM method for the equations (9.1.5) — (9.1.8). 
We are looking for two functions Ui 6 = 1,2 such that their weighted sum u = Xi u i + 


X 2 U 2 satisfies the Poisson equation (4.0.1) — (4.0.3) or its weak form 4.1.5 


This weak form of the iterative ASDDM now reads as: 


Algorithm 1 ASDDM Weak form 

Given u ° 6 = 1,2, such that u®\ r ; = g\r. t , i = 1,2, 

Define convergence level e 
while Error > e do 

Find u " E 'H 1 (Dj) : such that u™\ r j = u™~ 1 \ r j and uf |r\ = <7|r; 



gradv • gradufdfl 


Ifii 


vfdQ + / vhdT 


for any v E Ul(VLi),i,j = 1,2 ,i ^ y. 

end while 


11 Operator form of the algorithm 

In order to proceed to the finite dimensional space computations from this continuous settings we 
need to introduce a few operators. Assume that the solution uf on the D,; domain is sought in 
the finite dimensional space 5,. Given a function u” E our aim is to impose the Dirichlet 

boundary conditions on the D,; domain, which requires us to find a function g, E 5,; which will be 
close in some specific sense to the function Uj on the interface boundary P ; . 

We proceed as following: first we need to project the function u™ onto the boundary Tj to get 
a function g E £ 2 (Tl) representing the Dirichlet boundary conditions on this part of the boundary. 
We refer to this operator as a trace operator. 

The next steps depend on whether the space Si contains a function which is an extension of g or 
not. If it does we just take any extension of g satisfying the Dirichlet boundary conditions on the 
parts of the boundary T, and continue the solution process as in the regular IGA case, see section 
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(3.2.1). We refer to the operator building the extension as extension operator. If Si does not contain 
an extension of g we need to find a function in Sj that would approximate the g on the boundary 
Y\. We perform this approximation in two steps. First, we build an auxilary space 5m of functions 
that are just restrictions of the functions in Si onto the boundary T \. We then find among these 
functions the one gi E 5j|m that is closest to the g. Operator doing this is called approximation 

1 i 

operator. The process is then continued as before by applying the extension operator to gi to obtain 
9i £ Si. 

We now define the operators we have just introduced and rewrite the iterative algorithm in the 
operator notations. 

Loosely speaking the trace operator is defined as following: given a function w from a smooth 
enough functional space defined on Ylj the trace operator returns function defined on T) which 
coincides with w on Y) , i ^ j. 

In order to define the trace operator exactly we use the following 

Theorem 11.1 (Trace theorem) Let Yl C M n be bounded with piecewise smooth boundary. Then 
there exists a bounded linear mapping 

II : 7-{}(YL) -»• £ 2 (fl), (11.0.3) 

||n(u)||o, 9 n < C'||«||i,a, (11.0.4) 

such that n(it) = u\dn Vu E C 1 (fl). 

The proof can be found in H31 

Since we assume the domain boundaries to be piecewise smooth and the solutions belong to 
7-^ 1 (0y) spaces, the Trace theorem provides us with the definition of the trace operator: 


Pi 


Uj | r j■ 


Now assume that S t = 'H l {Yli) and the extension is possible. Given a function defined on T) we 
need to extend this function to the whole subdomain 11; while keeping the smoothness properties 
(almost everywhere with respect to the Lebesgue measure) of the function and in such a way that 
the new function will coincide with the Dirichlet boundary condition g on the outer boundary Tj. 
More precisely, the extension operators look like: 
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Vi ->■ Ui such that Ui\ T j = Vi,Ui\ Ti = g |r r 

Here we assume that the conditions Ui Lj = Vi,Ui |r !; = g |r, ; define some continuous function. 

This condition is satisfied throughout the whole discussion. 

Denote the corresponding billinear and linear forms on the subdomains as a,i,Li,i = 1,2. The 
continuous version of the ASDDM now reads as follows. 

Algorithm 2 ASDDM with Trace Operator 

Given u® G Ld(fh),* = 1,2, such that u® |r s = g\vi-,i = 1,2, 

Define convergence level e 
while Error > e do 

Find u ” G : suc/i that 

ai(ui - EiPiU™^ 1 ,Vi) = Li(vi ) - a^EiPiU^ 1 ,Vi) 


for any v t G Ul{VLi),i,j = 1,2,* ^ j. 

end while 


When we consider any iterative procedure the main question that arises is whether this process 
converges, and under which conditions. The convergence of this classical Schwarz algorithm was 
shown by P.-L. Lions m under some regularity assumptions. 

In the next section we consider the case when the approximation need to be performed before 
we can extend the boundary function to the whole domain and give a deeper overview of the trace 
operators we implemented. 

12 Discretized Schwarz Additive Algorithm and Trace Operators 

We now assume that we deal with the finite dimensional function spaces and S 2 instead of 
’H 1 (Di) and Ft 1 (^2), respectively. It means that we can not just project the function from Sj onto 
the boundary T 3 - and extend it to the whole fh domain, since this projection, generally speaking, 
will not belong to Si restricted to So we have to introduce an additional approximation operator 
as was explained in the previous chapter. This operator acts on functions from ’H 1 (T^) and maps 
them into the subspace SjLj CH^T-) of restrictions of functions from Si onto the boundary T|. 

More specifically, in our IGA environment we want this operator to build a B-spline or NURBS 
approximation of the function in hand. In the previous notations the approximation operator will 


26 





look as following: 


At : n\Ti)^Si\ rJ 

A i 

w —> v G <SLj, approximating w. 

L i 

where the approximating function can be thought of as minimizing distance to the function being 
approximated in some specific norm. 

With this notations the numerical algorithm will look like: 

Algorithm 3 Discrete ASDDM with Trace Operator 

Given initial guesses u® G <Sj, * = 1,2 such that tt? |r f = dhlr^ * = 1, 2 
Define convergence level e 
while Error > e do 
Find uf G S t such that : 

ai(u ” - EiAiPiU^ 1 ,Vi) = Lifa ) - ai{E i AiP i u 1 j~ 1 ,Vi) 


for any Vi G i,j = 1, 2, i ± j. 

end while 


Figure 12 illustrates the described ASDDM applied to a two domains problem. 


13 The Trace Operator 

In the present work we consider two ways of projection and compare their performance. Assume that 
both domains Pi and P 2 are B-spline surfaces with the corresponding piecewise invertible mappings 
Fi : Pi —> Pi and F 2 : P 2 —> P 21 where the parametric domains are Pi = P 2 = [0,1] x [0,1]- 


13.1 Exact trace operator 

We define the exact projetion operator as before: 

we get is a linear combination of 
the solution u'j 1 obtained after 


u j u j I r J ' 

% 

After each iteration on each domain the approximate solution 
images of B-spline (or NURBS) basis functions. Let us consider 
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Iteration n 


fll 


1 


Calculate the solution u x 



Project the solution u 2 onto the 
boundary of flj 


Approximate the boundary 
conditions by a NURBS curve 


Extend the NURBS curve to the whole 
domain fii 


New boundary conditions for fli 


Iteration n+1 


Calculate the solution u x 


Project the solution u 2 onto the 
boundary of fii 


T 


n 2 




Calculate the solution Ui 


Project the solution u 2 onto the 
boundary of fl 2 


Approximate the boundary 
conditions by a NURBS curve 


Extend the NURBS curve to the whole 
domain A2 


New boundary conditions for fl2 


Calculate the solution u- 



Project the solution u 2 onto the 
boundary of CI2 


T 


Figure 12: Operations performed on each iterative step in ASDDM with two subdomains Hi and H 2 
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Parametric spaces 



Figure 13: Domains and the pullback F ci ^ de ( r|) of the boundary r3> 


n — 1-th iteration on the subdomain Qj: 


u] 1 (x,y) = 4>k{x,y)uj1 \ (13.1.1) 

k&J 

where uj = {ujAiej is the vector of degrees of freedom. 

We want to evaluate this function at some point (x, y ) £ flj. Recall that the pull-back F ; - always 
exists since the geometrical mapping Fj that we use was chosen to be invertible: (£, y) = FT 1 (x,y). 
The problem here is that we cannot obtain a closed analytic form of FJ 1 , which means that in the 
real computations we will need to use some numerical methods to get the pre-image (£, rj) of the 
point (x,y). 

The resulting function is obtained by evaluating the solution w" -1 at any point (x,y) 

in this way. 


In the Figure 13.1 you can see an example of a pullback of the boundary T 2 into the parametric 
space Di 


13.2 Interpolation trace operator 

The second approach of performing the projection we used was the interpolation trace operator. 
Given the geometrical mapping F we can generate a mesh Ai in the physical space by applying 
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In-mesh interpolation 
(Physical space) 



this mapping to some mesh At in the parametric space. On the n-th iteration we can evaluate the 
solution u^j = 1, 2, in the physical domain f lj at the mesh points we have just constructed. Our 
aim is to approximate this solution at any point of the domain f lj. We obtain it by taking the values 
of the function u™ in the constructed mesh points and interpolating them into the points of interest, 
e.g. points on the curve T^. 

There is a wide choice of interpolation methods that we can use. In our research we used linear 
and cubic interpolations. The interpolation is obtained in the following way: we, first, divide the 
physical domain into triangles by applying the Delaunay triangulation algorithm to the mesh Ai. 
Then for any point ( x , y ) 6 Tl we define to which triangle it belongs and interpolate the value of 
the function u™ at this point using the values of this function in the vertices of the triangle we have 
chosen. 

Pi 

Uj ->• Iv{uj ), 

where Id is the interpolation operator, depending on the triangulation V of the mesh Ai. We will 
use the notation P l to denote the linear interpolation trace operator. 

In the next section we treat a one-dimensional example of ASDDM. We give an overview on the 
construction of the solver matrices and their structure. 
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Figure 15: One-dimentional Domain Decomposition 


14 One-dimensional ASDDM : Discretized Algorithm 


To illustrate the ASDDM on non-matching grids we consider here a simple one-dimensional Poisson 
equation on a line segment domain 0 = [0,1]. 

-u xx = f on Di, it(0) = u{ 1) = 0. (14.0.1) 


We decompose the domain into two overlapping domains Di = [0, (3) and D 2 = (a, 1], a < f3 such 


that Di U D 2 = D = [0,1] as shown in the Figure 14 We assume for convenience that the mapping 
F is the identity mapping, which implies that we may consider all the equations in the parametric 
spaces. We denote the knot vectors as Hi = [0,... /3] and H 2 = [cc,... 1] and the numbers of degrees 
of freedom as Nhi and , respectively. Also we assume that the knot vectors are uniform and open 
and denote the degrees of the polynomials p\ and P 2 on the subdomains Di and D 2 , respectively. To 
avoid too many indices we denote the solution on Di as v and the solution on D 2 as w. 

The two equations which we are going to discretize are: 


v xx = f onni,v n (0) =0,v n (/3) = w n 1 (/3). (14.0.2) 

w xx = / on D2,u; n (0) = 0 ,w n (a) = v n ~ l (a). (14.0.3) 


The numbering of the basis functions is exactly the same as we used in section [5] . 

Consider now the first subdomain D]. Only the first basis function cj>i is non-zero at 0 satisfying 
cf)\ (0) = 1. In the similar way at the second end of Di ipN hl (/3) = 1 and all the other basis functions 
vanish at point (3. It means, that in the one-dimensional case the imposition of Dirichlet boundary 
conditions becomes trivial: we just need to take the first and the last degrees of freedom equal to 
the projections of the solution on the other subdomain at the boundary points; the approximation 
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operators A % become identities. 

Denote the approximate solution on the first subdomain as v and its vector representation in the 
B-spline basis as v = Since the boundary condition at 0 is constantly zero, the first degree 

of freedom is 0: v\ = 0. The treatment of the boundary condition at the point (3 is more involved. 
We need to project the solution w n_1 of the previous iteration on the second subdomain 0,2 onto the 
boundary of fli which in our case consists of the only point (3. So we can write ujy = P\{w n ~ l ((3 )), 
where Pi is the trace operator. 

The same reasoning applies on 0,2- 

Denote the vector of degrees of freedom of the solution w as w. The solution w n ~ l now reads as 

N h 2 

i =1 

All the operators are acting on finite dimensional spaces, thus they have matrix representations in 
the basis we have chosen. Since the degree of B-splines on D 2 is P 2 we may conclude that there are 
no more than p 2 + 1 non-zero basis functions 1 pi at the point (3. As we discussed above, we may use 
different trace operators in order to project the solution on the subdomain O 2 onto the boundary 
T^ of the subdomain D]. In what follows we will consider both cases: the exact trace operator and 
the interpolation trace operator. 

We consider now the two kinds of trace operators we introduced in the previous chapter. 


14.1 Exact trace operator P e 


We take the value v 1 ^ of the boundary degree of freedom to be equal exactly to the value of the 
function w n ^ 1 (f3). For this trivial trace operator, the value we get is: 

I 0 \ 


v N hl = (0, ... ... ,ip i+P2 (l3),... ,0) ■ 


w 


n— 1 


W 


n— 1 
i +1 


W. 


,,n-l 

i+P2 


\ 0 ) 
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The operator matrix for the boundary of the fli subdomain is: 


^0 .. 

0 

0 

0 

. 0^ 

0 .. 

0 

0 

0 

. 0 

v° •• 

• MP) 

A+i(P) •• 

V’ i+p2 {P) 

• <7 


where Pf € W NhlxNh2 . 

Denote by A\ the stiffness matrix for the first subdomain. As before, we partition the sets of 
indices of basis functions J = { 1,2,..., N^j }, j = 1,2 into two subsets. The subset X C J of the 
inner degrees of freedom and B C J of the boundary degrees of freedom. In the one dimensional case 
B = { 1, Nhj} and X = dX\{ 1, Nhj } = {2,3,..., N^j — 1 },j = 1, 2. So the restriction of the stiffness 
matrix corresponding to the inner degrees of freedom reads as A\ (X, X). 


Since the basis functions are not local, when we impose the Dirichlet boundary conditions on 
some of the degrees of freedom, we need to subtract the function </>(£) from the solution v n . 
We do it by subtracting the corresponding values from the degrees of freedom of the basis functions 
k which domain intersects with the domain of (j>N hl - Now the discretized equation for the 
subdomain Di can be rewritten as: 


(\ 0 

0 o\ 


(o\ 


fl 

0^ 

0 Ai(I, X) 

0 

• v n = 

fl 

+ 


-Ai(l,B) 

\0 0 

0 l) 


\o) 


\0 

l) 


^0 ... 0 

0 ... 0 


0 ... 0 

0 ... 0 


o\ 


\0 ... A(P) A+i(P) • •• A+ P2 (P) ••• °/ 


(14.1.1) 


where the vector fi corresponds to the inner degrees of freedom of the first subdomain. Exactly the 
same reasoning can be applied to the second subdomain Q 2 to get the following discretized equation: 


(\ 0 

0 o\ 




fl 

0^ 

0 A 2 (X,X) 

0 

• w n = 


+ 


-A 2 (1,B) 

\0 0 

0 1 ) 


\0^ 


\0 

1 ) 


[0 ... (f>j(a ) (j) j+ i(a) ... (frj+vi ( a ) ••• 


: 


0 0 

0 0 


0 ... 0 

0 ... 0 J 


(14.1.2) 
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Let us denote the vector of degrees of freedom of the “whole” approximate solution u n as u n = 


Now, we can unite the two equations (14.1.1) and (14.1.2) into a single matrix equation 


w 


A • u n = f + A dir ■ P -u n \ 


(14.1.3) 


where 


A = 


(l ... 0 

0 A\{X.X) 0 
0 ... 1 


O 


\ 


with ii <E R NhlxNhl ,A 2 € 


Adir — 


(l ... 0 

—A\{I,B) 

0 


o 


\ 


o 

1 ... 0 
0 A 2 (Z,1) 0 
0 ... 1 ) 


/ 0 \ 
fl 
0 
0 
h 

w 


A 1 O 

o i 2 , 


/ = 


o 


1 


1 


0 


-A 2 {T,B) 

0 ... 1) 


■A-dirl O 

O -Adir2 . 


with Adiri G 'R NhlXNhl , Adir2 £ ]R Ar ^ 2XAr ^ 2 , 
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p = 


o 


0 

0 


0 

0 

A(P) 


0 

0 


V’j+P2 iP) 


cj)j(a) 


4>j+pi ( a ) 


o 



0 

0 


P 2 o 

with P{ G R NhlxNh2 ,P% e l^x Nh K 
The iterative scheme is given by: 



O Pf 


u n_1 = 


pe 

r 2 


O 


o\ 

0 


(14.1.4) 


(14.1.5) 


(14.1.6) 


We emphasize here that the matrix 

‘A x O 

O i 2y 

is not symmetric since the domains may have different number of elements, or even different orders 
of approximation. 


14.2 Linear interpolation trace operator 


We consider here the linear interpolation operator we discussed in the previous chapter. 

The algorithm works as was explained above: consider the first subdomain Oi = [0,/3). In order 
to construct a linear interpolation of the solution w n_1 at the point f] = f5 we need to find the knot 
span 1 ) which contains (5. We take the values of the function u ; n-1 at the ends of this chosen 

span interval uf n- 1 (? 7 *) and w 11-1 and their weighted sum: 


V N h 1 - 


Vi+i ~ Vi 


-w 


r.n—1 


(Vi) + 


Vi+1 - P , 
r)i +1 - Vi 


-w 


-:n— 1 


(Vi+ 1)- 
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Figure 16: One-dimensional Schwarz Method for the equation u xx = —1, using the exact trace 
operator P e and zero initial guesses 



Figure 17: One-dimensional Schwarz Method for the equation u xx = — 1, using the linear interpola¬ 
tion trace operator P l and zero initial guesses 
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We see that the linear interpolation trace operator is, actually, a convex sum of the two exact 
interpolation operators corresponding to the points rji and Vi+i'- 

P{ = P ~ Vi P e (Vi) + Vi+1 ~ f3 P e (Vi+i)- (14-2.1) 

v %+1 Vi n >+1 Vi 

Consequently, the matrix of this interpolation trace operator can be obtained as the convex sum of 
the matrices of the exact trace operators at the points rji and Vi+i- 


p i = P-Vi 

1 77j+i - rji 


+ 


Vi +1 ~ P 
Vi+I - Vi 


0 . 
0 . 

v° • 

o ... 
0 ... 


0 

0 


... 0 

... 0 


fpi(Vi) 4>i+l{vi) ... Tpi+itoiVi) ■■■ 0/ 

0 0 ... 0 ... 0 

0 0 ... 0 ... 0 


\o ••• A+i(vi+i) ^wivi+i ) iPi+P2+i(vi+i) ••• o 


The iterative scheme is now: 

' A x O 


n p i Adi r i O 

u n = f + 

O A 2 \ O Adir2 . 


O P\ 1 

1 1 • u ™- 1 = 

Pk O 


f + 


O Arlirl p} 


dir l r l | n—1 

li. , 


, AdiriPo O 


where 


0 +i 0 

and a belongs to the knot span [0,0+ 1 ) • 


“ ‘ ^ p'te) + 


pi = ,x ^ D e 


0 +i 0 


(14.2.2) 


(14.2.3) 


(14.2.4) 


14.3 Initial guesses 


Of course, there is some freedom in choosing the initial boundary values u°(/3) and w°(a). It turns 
out that the method converges to the same approximate solutions regardless of the initial guesses 
since the convergence is guaranteed by the properties of the matrix, which is a M-matrix in our case 
as shown in the Appendix. 


You can see an example of a convergence procedure for non-zero initial guesses in the Figure 18 
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Figure 18: Convergence of the one-dimensional SADDM with non-zero initial guesses. Here the 
initial values for the inner boundary points are taken v°(/3) = /3 and -u; 0 (a) = a 

14.4 Multi-dimensional domains 

When we consider a multi - dimensional case we have to perform an approximation of the boundary 
conditions since the boundary condition function is be a B-spline or NURBS curve or surface in 
general. This approximation can be performed in one of the ways we have already discussed above 
in section [7] : quasi-interpolation or least-squares approximation. All the other steps of the solution 
are similar to the one-dimensional case. 

In the next chapter we consider two and three dimensional applications of ASDDM as well as its 
parallelisation. 

15 Numerical Results 

In this section we provide illustration of the convergence properties of the ASDDM in our imple¬ 
mentation. As we have already mentioned one of the main advantages of the ASDDM methods is 
that the meshes of different patches are non-matching. This important feature makes it possible to 
apply this method in CSG and zooming problems. In particular, it allows to perform local zooming 
at places where the solution possesses weak singularities. Below we demonstrate some examples of 
this technique and its usage. We also describe how we implement the parallelized version of the 
solver and show a few three-dimensional multi-patched examples obtained with it. 
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Figure 19: £ 2 -error iterative convergence depending on the polynomial degree, one-dimensional case. 


As already mentioned before in all the model problems considered in this work the errors are 
£ 2 -norm differences between the exact solutions and the approximations, see section (4.1.2). 


16 One-dimensional examples and Analysis 


16.1 Dependence of the convergence on the B-spline degree 


First, we treat the one-dimensional problem to show how the iterative process depends on the 
degree of the B-splines used. The graph given in Figure 16.1 shows the £ 2 -errors as a function of 


the iteration number for different polynomial degrees. Obviously, higher polynomial degrees lead to 
better precision of the final result, but the rates of convergence remain the same which makes the 
iterative procedure time consuming. 


16.2 Dependence of the convergence on the overlapping area 

Many times when the Domain Decomposition technique is applied to the original domain we have 
some freedom in choosing the subdomains. For example, when a local zooming is considered the 
only thing which is predefined is the domain we zoom in, but the choice of the second subdomain is 
given to us. The question is how the relative overlapping area influences the convergence process. 
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“overlap 0.5 
“overlap 0.4 
overlap 0.3 
-overlap 0.2 
overlap 0.1 


L2 error depending on the overlap area 


0.0001 


0.00001 


iterations 


Figure 20: £ 2 -error iterative convergence depending on the ovrelapping area 


In Figure 20 we compare different domain decompositions for the same one-dimensional example. 
The overlapping is measured as the ratio of the overlapping area to the area of one subdomain. 
Apparently, the larger the overlapping is, the faster the iterative convergence becomes. An open 
question is how would one balance the overlapping and the number of degrees of freedom versus the 
number of iterations to optimize the computation. 


16.3 Dependence of the approximate solution on the mesh size 

Figure [21] shows how the approximate solution depends on the refinement for different degrees of 
B-splines. You can see that the dependence is close to theoretically predicted, given by the formula 

m- 


inf ||u o F — s ||£2 < Ch p+1 \u o F|-^ P +i, (16.3.1) 

seV 

where u is the exact solution and C is a constant that may depend on the degree p of the polynomials. 
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Figure 21: One-dimensional £ 2 -error of the approximate solution depending on the mesh size for 
different polynomial degrees 


17 Convergence properties in multi-dimensional cases: Applica¬ 
tion to Zooming 


The formula 16.3.1 becomes different in the multi-dimensional case, but still the approximated 
solution depends on the mesh size in almost the same way m- 

Consider the following two-dimensional problem: 

Let P be an open circle of radius 3 placed at the origin: 


— A u = sin(x 2 + y 2 — 9) — 4cos(x 2 + y 2 — 9) on P, 


(17.0.2) 


"w|ar2=o- 


(17.0.3) 

The exact solution is u = sin(x 2 + y 2 — 9) 

Apply zooming by considering P as a union of an annulus and a square. P = P annulus U P square 


as shown in the Figure 22 


In Figures 23 and |24| below you can see the convergence properties and the numerical solutions 
we obtained in this case by applying the ASDDM. 
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Figure 23: Annulus £ 2 -errors depending on the mesh size for different polynomial degrees 
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Figure 24: Numerical solution 


18 Two-dimensional examples 

18.1 Non-matching domain meshes 

An example of non-matching meshes and a solution obtained by the Domain Decomposition method 
is shown in Figure[25j The domain is taken from the logo of the Domain Decomposition organization, 
seeFigure m We also give in Table 1 the reaction between overlap and DD iterations. 


18.2 Application of zooming: weak corner singularity 


We now consider some engineering applications of the Domain Decomposition techniques used for 
zooming. One of the main applications of ASDDM are the problems where the solution possesses a 
weak singularity at the corner. The example we treat here is taken from [7]. 

Given the following domain: 


P = (0, 0) € M|0 < p < 3; < 0 < it}, 


(18.2.1) 


which is an open circular sector of angle a 


^ as shown in Figure 


26 


and the 


equation: 


A u = 0 in P, 


(18.2.2) 


0 : on Tj cTij = <9P 

0(a-0) :onr 2 = r n \r 1 . 
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Figure 25: Numerical solution over the DD logo. 



Figure 26: Problem with a weak singularity at the corner. The domain. 
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Table 1: Number of Iteration and max error by Overlap Distance 


Overlapp 

Iterations 

0.1 

max error 

0.25 

at each iteration 

0.5 

0.75 

1.0 

1 

.6 * 10 _1 

.9* 10' 1 

1.2 * 10 -1 

1.2 * 10 -1 

0.7* 10" 1 

2 

.3 *10 -1 

.3* 10” 1 

0.16 *10 -1 

0.12* 10” 1 

0.11 *10” 1 

3 

.15* 10" 1 

.024 * 10" 1 

0.01 * 10" 1 

0.01 * 10” 1 

0.022* 10” 1 

4 

.05 *10 _1 

.08* 10" 1 

0.034* 10” 1 

0.005 * 10” 1 

0.006 * 10 -1 

5 

.03 * 10" 1 

.007* 10" 1 

0.006* 10” 1 

0.006 * 10" 1 

0.001 * 10" 1 

6 

.02 *10 -1 

.004 * 10 -1 

0.002 * 10” 1 

0.003 * 10 -1 

c 

7 

.02 * 10" 1 

.002 * 10" 1 

0.0015 * 10" 1 

0.002 * 10" 1 

c 

8 

.01 *10 _1 

.002 * 10 _1 

0.001 * 10" 1 

0.001 * 10 _1 

c 

9 

.006 * KU 1 

.001 *10" 1 

c 

c 

c 

10 

.003 * IQ -1 

c 

c 

c 

c 


The exact solution is given by the following uniformly convergent series: 


Uex(p,8) = - 4 

7r ' n 6 

n= 1 , 3 , 5 ,... 


2 n 
3 


Sill 


(?) 


(18.2.3) 


We see from (18.2.3) that the radial derivative is unbounded when approaching the origin which 
may cause the numerical solution to be unstable at the vicinity of the origin. In order to zoom the 
solution near the corner we decompose the domain in the following way: D = D 0 U jUD 200m , as shown 
in Figure [27| 

Figure [28] shows the numerical solution and the radial derivative calculated for this problem. The 
principal order of the radial derivative of the exact solution is — and the numerically computed 



Figure 27: NURBS definition, Boundary and Domain Decomposition 
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Figure 28: Numerical solution of (3.2) and the corresponding radial derivative 


result using the Domain Decomposition technique coincides with this value. 

We have tested different forms of the zooming region as shown in Figure [29] and in all the cases 
the results coincided. 

The solution in this case is not regular and possesses a weak singularity at the corner and the 


error estimate 16.3.1 does not apply anymore. The dependence of the approximate solution on the 
mesh size becomes different from the regular case and the rate of convergence does not improve with 


greater polynomial degree as shown in Figure 30 


19 Three-dimensional examples 


In this section we present some examples obtained for multi-patched bodies. Consider the equation: 


— A u = sin(x + y + z) on D, 


(19.0.4) 


u\da. = sin (x + y + z), 


(19.0.5) 


where 9 is a chain of 5 overlapping cubes, as shown in Figure 31 The exact solution is u = 


sin(x + y + z). You can see the rates of convergence of the solution for all the 5 cubes in the Figure 
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“degree 2 
“degree 3 
-degree 4 


Figure 30: £ 2 -error of the approximate solution depending on the mesh size for different polynomial 
degrees 
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Figure 31: Chain of cubes with the exact solution sin (a: + y + z) 



Figure 32: Iterative £ 2 -error for the chain of cubes 
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OPEN POOL 



Synchronization 


Figure 33: Synchronization algorithm 

19.1 Parallelized algorithm and multi-patched problems 

As we have mentioned in the beginning of this work, one of the main advantages of the ASDDM is 
that it may be easily parallelized. The modern computers are all multi-core which provides us with 
an opportunity of real parallelization of the code, meaning that the solution on each patch may be 
computed independently and simultaneously on a different processor. 

The algorithm we implemented works as following. Before we start the iterations we precompute 
all the date that the solver needs, then we open a worker thread for each subprocess, corresponding 
to a patch, and start the iterative procedure. On each iteration the n-th thread solves the equation 
on the n-th subdomain according to the scheme shown in Figure [12] of the previous chapter. After 
the n-th iteration computation is finished, the data is synchronized between the threads. This 
synchronization is necessary to impose the correct boundary conditions for each subdomain on the 
next iteration. 

The schematic synchronization mechanism is shown in Figure [33| 


49 














Figure 34: Elasticity example 

19.2 Elasticity example for multi-patched domain 

As an example of parallel computation we took an elasticity problem. The domain is a hollow thick 
half-ring with fixed ends which is placed into a gravitational field. This patch was decomposed 
into 8 overlapping subdomains as shown in the Figure [34] and the solution was obtained using the 
parallelized method described above. 

We performed numerical tests on the 4-processors machine and achieved a full 4-times enhance¬ 
ment in the time of computations. 

19.3 Statistics 

In this section we illustrate the dependence of the number of degrees of freedom in the problem as a 
function of dimension, refinement and polynomial degrees. Of course, these numbers do not depend 
on the isogeometrical mapping, so we calculated the data using a cube domain. 

In the two-dimensional case we obtained: 
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degree of the B-splines 


Partitions of the mesh 2 3 4 


0 

10 

27 

52 

1 

27 

52 

85 

2 

52 

85 

126 

5 

175 

232 

297 

10 

540 

637 

742 

11 

637 

742 

855 

20 

1870 

2047 

2232 


In the three-dimensional case we obtained: 


degree of the B-splines 
Partitions of the mesh 1 2 

1 765 624 

2 2100 1275 


The number of degrees of freedom, and consequently, the computational time grows rapidly 
with the dimension of the problem, especially if the mesh need to be refined. Here the Domain 
Decomposition techniques become very helpful. 

20 Conclusion 

The main purpose of this work was to establish the connection and application of the recently 
discovered IGA framework to solve PDEs using Additive Domain Decomposition Method. It turns 
out that IGA may be successfully used for solving different problems even on complex geometries, 
thus giving us all the advantages of IGA without restrictions on the geometry. Another important 
application of these methods are zooming problems. The IGA framework is restricted to its native 
techniques like T-splines and here the Domain Decomposition becomes very useful. 

More research is needed is to explore ways of approximating and imposing the Dirichlet boundary 
conditions. Different techniques of imposing boundary conditions may influence the computational 
times and the precision of the approximate solution significantly hence a research subject is improving 
techniques to define such BCs on trimmed patches. 
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The question of optimal overlap choice ( a large overlap seems better, but more costly in degrees 
of freedom) needs more study. Also what is the best way to define the solution on the intersections 
is an opened question. 

Another important direction for research is the preconditioning of the solver in large scale prob¬ 
lems. An efficient pre-conditioner may improve the computational times and simplify the structure 
of the solver( and insure convergence , since we do not have a maximum principle here.) 

Finally a CSG construct is not made only of unions of primitives, and it my be the results 
of Boolean differences, substractions etc.,thus DD should be coupled with other methods such as 
fictitious domains . 
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